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We present numerical simulations based on a Boltzmann kinetic model with compet- 
ing interactions, aimed at characterizating the rheological properties of soft-glassy 
materials. The lattice kinetic model is shown to reproduce typical signatures of 
driven soft-glassy flows in confined geometries, such as Herschel-Bulkley rheology, 
shear-banding and histeresys. This lends further credit to the present lattice ki- 
netic model as a valuable tool for the theoretical/computational investigation of 
the rheology of driven soft-glassy materials under confinement. 
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1. Introduction 

Many applications in modern science, engineering and biology have prompted the 
recent blossoming of research in the rheology of soft-flowing and noncrgodic materi- 
als, such as emulsions, foams, gels (Larson 1999,Coussot 2005, Chaikin & Lubensky 
1995, Lyklema 1991, Evans & Wennerstrm 1999, Degennes 1979, Doi & Edwards 
1986, Grosberg & Khokhlov 1994, Weaire & Hutzler 1999). The theoretical un- 
derstanding of such materials raises interesting questions on its own. Indeed, since 
soft-materials share simultaneously many distinctive features of the three basic 
states of matter (solid, liquid and gas), their quantitative description does not fall 
within the traditional methods of equilibrium and/or non-equilbrium statistical 
mechanics. New concepts and theoretical paradigms are required to adjust many 
non-standard features, such as rheology, long-time relaxation, dynamic disorder, 
ageing and related phenomena. In particular, under simple shear conditions, some 
of these complex fluids may separate into bands of widely different viscosities. This 
phenomenon, known as "shear banding" (Berret 2005), involves inhomogeneous 
flows where macroscopic bands with different shear rates or shear stresses coexist 
in the sample. Although shear banding is attributed to a shear-induced transition 
from a microscopic organization of the fluid structure to another, it still raises lots 
of theoretical and experimental challenges (Manneville et al. 2007, Sollich et al. 
2009, Fardin et al. 2010). 

Simultaneously, the need for new tools of analysis also emerges for computational 
studies, which typically involve many interacting space and time scales, the latter 
case being particularly acute in view of the aforementioned long-time relaxation. In 
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the recent past, we presented a new conceptual/computational scheme for the nu- 
merical simulation of soft-flowing materials (Benzi et al. 2009, Bernaschi et al. 2009, 
Benzi et al. 2010). The scheme is based on a (Lattice) Boltzmann (LB) formulation 
(Bathnagar, Gross & Krook 1954, Benzi et al. 1992, Chen & Doolen 1998, Gladrow 
2000) for interacting binary fluids (Shan & Chen 1993,1994), in which, by a proper 
combination of short-range attraction and mid-range repulsion (competing self in- 
teractions) , an effective form of frustration was encoded within the physics of the 
binary mixture (Benzi et al. 2009). More specifically, by tuning the above interac- 
tions in such a way as to bring the surface tension down to nearly vanishing values, 
many typical signatures of soft-glassy behaviour, such as long-time relaxation, dy- 
namical arrest, ageing, and non-linear Herschel-Bulklcy rhcology, have been clearly 
detected in numerical studies (Benzi et al. 2009, Bernaschi et al. 2009, Benzi et al. 
2010). In this paper, we further elaborate along these lines by performing numer- 
ical simulations of confined flows and looking at their rheological properties, i.e. 
the relation between the applied shear and the developed stress by the fluid. By 
performing numerical simulations over a wide range of shear rates, we detect the 
emergence of shear banding effects, i.e. a stress plateau, which is shown to corre- 
spond to a wide range of stationary shear rates. These results lend further credit 
to the Boltzmann kinetic formulation as a valuable tool for the systematic study of 
the emergence of non linear rheological properties from mesoscopic interactions. 



2. The Model: Lattice Boltzmann equation with competing 

self interactions 

In this section we review the main properties of the lattice Boltzmann model used in 
the numerical simulations. Further details can be found in our recent work (Benzi et 
al. 2009). The starting point is the lattice transcription of the generalized Boltzmann 
equation (Bathnagar, Gross & Krook 1954, Benzi et al. 1992, Chen & Doolen 1998, 
Gladrow 2000) for a multicomponent fluid with S species, as inspired by the work 
of Shan & Chen (Shan & Chen 1993, Shan & Doolen 1995): 

f is (r + CiAt,t + At) - f is (f,t) = - — [f is (f,t) - f^ q \p s ,u + T s F s /p s )}. (2.1) 

Ts 

In the above, fi S (f , t) is the probability density function of finding a particle of 
species s — 1...S at site r and time t, moving along the i-th lattice direction, 
defined by the discrete speeds c\ with i = 0...b (sec Figure 1). For simplicity, the 
characteristic propagation time lapse At is taken equal to unity (At = 1) in the 
following. The left hand-side of (|2.1[) stands for molecular free-streaming, whereas 
the right-hand side represents collisions as a simple time relaxation towards the 
local Maxwellian equilibrium f^ q \p s ,u) on a time scale r s . The local Maxwellian 
is truncated at second order, an approximation which suffices to recover the correct 
athcrmal hydrodynamic balance 

ft\ P s,u) = ^Ps (l+ ^ + {C - C X c f s5ab) U aUb ) 

with c| the square of the sound speed in the model and 5 a b the Kronecker delta 
with a, b indicating the Cartesian components (repeated indices are summed upon). 
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The w^ eq,> 's are equilibrium weights used to enforce isotropy of the hydrodynamic 
equations (Benzi et al. 1992, Gladrow 2000, Chen & Doolen 1998). The equilibrium 
for the s species is a function of the local species density 

Ps(r,t) = ^2f is (r,t) 

i 

and the common velocity, defined as 

The common velocity receives a shift from the force F s acting on the s species (Shan 
& Chen 1993, Shan & Doolen 1995). This force may be an external one or due to 
intermolecular interactions. The pseudo-potential forces embed the essential fea- 
tures to achieve phase separation (non-ideal equation of state and non-zero surface 
tension), as well as a mechanism of frustration through competing self- interactions. 
More specifically, within each species, the forces consist of an attractive (denoted 
with (att)) component, acting only on the first Brillouin region (belt, for simplic- 
ity), and a repulsive (denoted with (rep)) one acting on both belts, whereas the 
force between species (X) is short-ranged and repulsive. In equations: 



where 



F s (i>,t) = Fj att) (r , t) + i?<™*> (f, t ) + (f , t) 
F( att \r 1 t) = -G ( s att ^ s (r,t) £ ( 2 - 2 ) 



iebeltl 



F^(r,t) = -G^ p Hs(r,t) £ p^ s {r h tfc - G< pe *ty.(f, t) ]T p^ s {n,t)c\ 

iebeltl iebelt2 

(2.3) 

^(r,t) = -—^-p s (f,t)J2 E G ss' Wip s >{fi,t)ci. (2.4) 

KPO ) s'^s iebeltl 

In the above, the groups "belt 1" and "belt 2" refer to the first and second Brillouin 
zones in the lattice and Cj, Pi,Wi are the corresponding discrete speeds and asso- 
ciated weights (sec figure [1] and table [1]). The interaction parameter G ss > — G S ' S , 
s' ^ s, is the cross-coupling between species, po is a reference density to be defined 
shortly and, finally, fi = 7 r + c\ are the displacements along the i-th directions. 
These interactions are sketched in Figure [1] for the case of a two component fluid 
(say species A and B). This model bears similarities to the next-to-nearest-neighbor 
frustrated lattice spin models (Shore & Sethna 1991, Shore, Holzer & Sethna 1992). 
However, in our case, a high-lattice connectivity (Sbragaglia et al. 2007) is required 
to ensure compliance with macroscopic non-ideal hydrodynamics, which is at the 
core of the complex rheology to be discussed in this work. To this purpose, the first 
belt is discretized with 8 speeds, whereas the second with 16, for a total of b = 25 
connections (including rest particles). All the weights take the values illustrated in 
Table □ 

The present model has already proven capable of reproducing several signatures 
of soft-glassy behaviour, such as long-time relaxation, ageing, and non-Newtonian 
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1 st belt 2 nd belt 
22 15 10 14 21 




23 18 12 19 24 



Figure 1. The discrete 25-speed lattice. Both belts are illustrated with the corresponding 
velocities. Further details can be found in the paper by Benzi and coworkers (Benzi et al. 
2009). 



Forcing Weigths (for Fj' ep ^) 



Pi 


= 247/420 


i = 


pi 


= 4/63 


i = l,4 


Pi 


= 4/135 


i = 5,8 


Pi 


= 1/180 


i = 9,12 


Pi 


= 2/945 


i = 13,20 


Pi 


= 1/15120 


i = 21,24 



Forcing Weights (for F { s att) and i?f ) 


w, = 4/9 


i = 


w s = 1/9 


i = 1,4 


= 1/36 


i = 5,8 



Table 1. Links and weights of the two belts, 25-speeds lattice. 



rheology. As a general remark, the onset of non-trivial behaviour has been found 
to associate with the regime of very-low surface tension, a physical quantity that 
the present scheme allows to tune through the coupling parameters G's and po. In 
particular, once the couplings G's are fixed, surface tension can be brought down to 

nearly vanishing values by increasing the value of po, because the repulsive intra- 

-. x 

species force, F s , contributing a positive surface tension, are scaled by a factor 
1/Pq, see eq. (|2.4[) . Full details can be found in our previous publications. To date, 
the aforementioned effects have been explored only in homogeneous, boundary-free 
geometries. Given the importance of relating to experimental results, it is of great 
interest to assess whether the above phenomenology can also be reproduced for the 
case of soft-glassy flows under geometrical confinement. This is precisely the task 
undertaken in the following section. 
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Figure 2. Stress vs shear for the case po = 0.70. The solid line is a linear fit (exponential 

in log representation). 



We have simulated a two-dimensional channel flow driven by the upper wall, mov- 
ing at speed U w along the mainstream direction, x. The channel measures L = 256 
lattice sites in length and H — 512 lattice sites in width. The initial conditions corre- 
spond to a modulated density, pa(x, y,t = 0) = p m (l + 0.1 sin (6.28^) sin (6.28^)), 
p B (x,y,t = 0) = p m (1 + 0.01 sin (6.28^) sin (6.28^)) with p m = 0.612 and zero 
flow. Periodic boundaries are imposed at the inlet and outlet sections. The coupling 



parameters are as follows G ( ° U) = -9.0, G ( ™ p) = 8.1, G ( ° tt) = -8.0, G% ep) = 7.1, 



and Gab = 0.405. The time relaxation is r = 1 for both species, corresponding 
to a kinematic viscosity v = 1/6 in lattice units. Under ordinary (Newtonian) flow 
conditions, the wall drive generates a linear flow profile u x (y) = U w y/H, with an 
associated homogeneous stress a xy = pU w /H, p = pv being the dynamic viscosity 
of the fluid. This is precisely the situation our model is expected to reproduce in 
the regime of sufficiently high-surface tension. 



In Figure [51 we report the average stress cr xy , as collected over the set of slices at 
x = 1, L, as a function of the applied wall-shear S w = U w /H. The simulations have 
been performed with po = 0.70, corresponding to a sizeable value of surface tension 
(7 ~ 0.05), for which Newtonian behaviour is expected. As one can appreciate, the 
stress-strain relation, in log-scale, is well fitted by an exponential curve, indicating 
a linear shear-stress relation. 

(b) Low surface tension: non-newtonian rheology and shear banding 

In order to investigate the emergence of non-Newtonian behaviour, we have 
repeated the simulations with a higher p = 0.83, corresponding to a much lower 
surface tension 7 ~ 5 x 10~ 3 . The corresponding stress-shear relation is shown in 
figurc[3] Several comments are in order. At low shear, S w < 2 10 -6 , the stress-shear 
still follows a linear, newtonian relation, as in the high-surface tension case described 



3. Numerical Results 



(a) High surface tension: Newtonian rheology 
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Figure 3. Left: Stress vs shear for the case po = 0.83. As we observe, there is a range 
of shear rates at which a Plateau in the stress tensor is observed. All the results have 
been obtained by averaging over many configurations and starting from the same initial 
condition. Right: we repeat the same measurements for the Stress vs shear starting from 
an initial configuration at high shear rates and decreasing the shear at the wall. Hysteresis 
effects are visible. 



previously. Upon increasing the driving shear, however, no further increase of the 
stress is observed, up to S w = 2 1CP 5 . This is the so-called shear-banding effect, 
signaling qualitative rearrangements in the structural response of the flow pattern, 
to be discussed shortly. By further increasing the driving shear, the stress curve 
regains a linear (exponential in log-scale) behaviour, corresponding however to a 
smaller, value of the effective viscosity. The emerging picture is that of a yield- 
stress, shear-thinning fluid, whose effective viscosity can be mapped into a Herschel- 
Bulkley relation of the form a xy = A + B S%, with A = 1 .48 x 1(T 6 , B = 0.000441 
and a ~ 0.25. The appearance of a non-zero yield stress, cry = A, prepares the 
ground for the emergence of shear banding. Indeed, the fluid can only flow in regions 
where the stress exceeds uy (the shear bands), while in the rest of the domain the 
fluid stands basically still, like a solid. 

(c) Low surface tension: hysteresis 

We have also investigated the occurrence of memory- dependent phenomena, such 
as hysteresis. To this purpose, we have repeated a series of simulations taking as 
an initial condition the steady state at S w = 9.76 x 10~ 5 , and then scanning S w 
backwards to decreasing values. As again apparent from figureEl the corresponding 
values of the stress do not lie on the same curve as before, but trace in fact a lower- 
lying curve. A similar behaviour is observed by repeating the same backward scan, 
starting from a higher value S w = 1.5 x 10~ 4 , although no crossover is observed in 
this latter case. By and large, all of the non-Newtonian features portrayed in this 
figure reveal a striking similarity with recent experiments on soft-glassy flows in 
microchannels (Manneville et al. 2007). 

(d) Spatial velocity profiles 

Finally, we inspect the spatial distribution of the average flow velocity, u x (y), 
inside the channel, see figure |U From this figure, it is apparent that in the low- 
surface tension, non-Newtonian regime, the velocity profiles tend to form a steep 
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Figure 4. Velocity profiles for po = 0.83 and different wall velocities corresponding to the 
following shears S w = 9.76 x 1(T 6 , 1.46 x 1CT 5 , 1.95 x 10" 5 . The vertical coordinate has 
been normalized between and 1 and the flow velocity has been normalized with respect 
to the corresponding wall velocity. 

layer next to the two walls (the shear bands), with a susbtantial flattening in the 
bulk region of the flow, where little dissipation takes place. This flattening is re- 
sponsible for the plateau in the shear-stress relation. Much like a turbulent flow, 
under increasing shear, the non-newtonian fluid reacts in such a way as to confine 
most of the shear, and stress, to an increasingly thinner near- wall layer. This is 
again in line with recent experimental observations (Manneville et al. 2007). 

4. Conclusions 

Summarizing, we have shown that the two-component Lattice Boltzmann model 
with competing interactions is capable of reproducing distinctive non-Newtonian 
features of driven soft-glassy flows in confined geometries, such non-Newtonian rhe- 
ology, shear-banding and hysteresis. The simulation data bear striking similarities 
with recent experimental results of driven soft-glassy flows in Couette geometries. 
Future work shall be directed to the quantitative comparison with experimental 
data, as well as to the investigation of the implications of low surface tension on 
the global and local rheology of soft-glassy materials. 
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